
%% plot data and secular variation

S = dir('/Users/manojnair/data/obs_mag_data/BGS data/*.nc');


selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
    'BLC';	'BMT';	'BOX';	'BOU'; 'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
    'EBR';	'EYR';	'FCC';	'FUR';	'GNA'; 'GLN';	'HLP';	'HRN';	'KAK';	...
    'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
    'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
    'STP';	'STJ';	'THY'; 'TUC';	'UPS';	'VIC'];

time_array_hour = (datenum(2010,1,1,0,0,30): (1/(24)): datenum(2014,12,31,23,59,30))';

for i = 1:length(S),
    filestring(i,:) = S(i).name(1:3);
end;

[data_obs, ia, ib ] = unique(filestring,'rows');

[spline_obs, ia, ib]= intersect(data_obs, selected_obs, 'rows');

% plot

f = figure(1);
set(f,'Position',[40          19        1557         937]);


annual_time_1 = datenum(2010,1:48,15);
annual_time_2 = datenum(2011,1:48,15);
actual_sec_time = datenum(2010,7:54,17);
annual_time = datenum(2010,1:60,15);

for i = 1:length(spline_obs),
    
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS data/' spline_obs(i,:) '_1995_2014.nc'];
    
    
    display(ncfname);
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, 131496, 43824, 0);
    
    if any(x_data),
        
        
        
        % fit splines
        L = isnan(x_data);
        
        data_array = x_data(~L);
       first_index = find(~L, 1, 'first');

        last_index = find(~L, 1, 'last');
        
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp_x=spline(b1,x_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_y=spline(b1,y_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_z=spline(b1,z_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            
            if  obj.geospatial_lon > 180,
                obj.geospatial_lon = obj.geospatial_lon - 360;
            end;
            
              
            icount = 1;
            for decyear = 2010:0.1:2015.0,
                [mag_igrf11(icount,:)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
                
                geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
                geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km                
                [mag_igrf12(icount,:)] = magsynth(geoc_r,geoc_lat,obj.geospatial_lon,decyear);
                icount = icount + 1;
            end;
            
            % secular variation is linear and same through out the model so
            % no need to calculate it through the loop.         

            
                geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
                geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km                

            
            sec_igrf12 = magsynth(geoc_r,geoc_lat, ... 
                obj.geospatial_lon,2015.5)-magsynth(geoc_r,geoc_lat, ... 
                obj.geospatial_lon,2014.5);

            sec_igrf11 = (igrf11magm(0,obj.geospatial_lat, ... 
                obj.geospatial_lon,2015.0)-igrf11magm(0, ... 
                obj.geospatial_lat,obj.geospatial_lon,2010))/5;
            
            sec_x_obs =  ppval(sp_x,annual_time_2) - ppval(sp_x,annual_time_1) ;
            sec_y_obs =  ppval(sp_y,annual_time_2) - ppval(sp_y,annual_time_1) ;
            sec_z_obs =  ppval(sp_z,annual_time_2) - ppval(sp_z,annual_time_1);
            
            
            
            
            % remove non data region.
            
            sec_x_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
            sec_y_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
            sec_z_obs(actual_sec_time > time_array_hour(last_index)) = NaN;

            
        
            
            f = figure(1);
            set(f,'Position',[ 560    39   778   909]);
            subplot(321);
            
            plot(actual_sec_time, sec_x_obs,'LineWidth',2);
            set(gca,'FontSize',16);
            hold on;
            plot(datenum(2010,1,1), sec_igrf11(1),'r*','MarkerSIze',20);
            plot(datenum(2014,4,20), sec_igrf12(1),'k*','MarkerSIze',20);
            h=legend('OBS','IGRF11','IGRF12');
            set(h,'FontSize',8,'Location','NorthWest');
            
            
            datetick('x','keeplimits');
            ylabel('X(nT/yr)');
            title(sprintf('%s lat %5.2f lon %5.2f ',spline_obs(i,:), ... 
                obj.geospatial_lat, obj.geospatial_lon));
            
            subplot(323);
            plot(actual_sec_time, sec_y_obs,'LineWidth',2);
            set(gca,'FontSize',16);
            hold on;
            plot(datenum(2010,1,1), sec_igrf11(2),'r*','MarkerSIze',20);
            plot(datenum(2014,4,20), sec_igrf12(2),'k*','MarkerSIze',20);
            
            
            datetick('x','keeplimits');
            ylabel('Y(nT/yr)');
            
            subplot(325);
            plot(actual_sec_time, sec_z_obs,'LineWidth',2);
            set(gca,'FontSize',16);
            hold on;
            plot(datenum(2010,1,1), sec_igrf11(3),'r*','MarkerSIze',20);
            plot(datenum(2014,4,20), sec_igrf12(3),'k*','MarkerSIze',20);
            
            
            datetick('x','keeplimits');
            ylabel('Z(nT/yr)');
            
            
            
             subplot(322);
            
            plot(time_array_hour,x_data);
            hold on;
            spline_fit = ppval(time_array_hour,sp_x);
            spline_fit(1:first_index) = NaN;
            spline_fit(last_index:end) = NaN;
            plot(time_array_hour,spline_fit,'r');
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf11(:,1),'r.-','MarkerSize',15);
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf12(:,1),'b.-','MarkerSize',15);
            
            h=legend('OBS','FIT','IGRF11','IGRF12');
            set(h,'FontSize',8);
            set(gca,'FontSize',16);
            
            ylabel('X nT');
            axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
            datetick('x','keeplimits');
            
            
            subplot(324);
            
            plot(time_array_hour,y_data);
            hold on;
             spline_fit = ppval(time_array_hour,sp_y);
            spline_fit(1:first_index) = NaN;
            spline_fit(last_index:end) = NaN;
            plot(time_array_hour,spline_fit,'r');
            
            hold on;
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf11(:,2),'r.-','MarkerSize',20);
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf12(:,2),'b.-','MarkerSize',20);
            
            
            set(gca,'FontSize',16);
            ylabel('Y nT');
            axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
            datetick('x','keeplimits');
            subplot(326);
            
            plot(time_array_hour,z_data);
            hold on;
            spline_fit = ppval(time_array_hour,sp_z);
            spline_fit(1:first_index) = NaN;
            spline_fit(last_index:end) = NaN;
            plot(time_array_hour,spline_fit,'r');
            
            hold on;
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf11(:,3),'r.-','MarkerSize',20);
            plot(datenum(2010:0.1:2015.0,0,0), mag_igrf12(:,3),'b.-','MarkerSize',20);
            
            
            set(gca,'FontSize',16);
            ylabel('Z nT');
            axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
            datetick('x','keeplimits');
            
            
            set(f,'PaperPositionMode','auto');
            saveas(f,['/Users/manojnair/projects/igrf12_validation/' sprintf('%s',spline_obs(i,:)) '_Sec_Var'],'epsc');
            
            
           % pause;
            
            clear  sp_x sp_y sp_z;
            clf;
            
        end;
        
    end;
    
end;


